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Calculating an accurate nutation time constant for a spinning spacecraft is an important 
step for ensuring mission success. Spacecraft nutation is caused by energy dissipation about 
the spin axis. Propellant slosh in the spacecraft fuel tanks is the primary source for this 
dissipation and can be simulated using a forced motion spin table. Mechanical analogs, such 
as pendulums and rotors, are typically used to simulate propellant slosh. A strong desire 
exists for an automated method to determine these analog parameters. The method 
presented accomplishes this task by using a MATLAB Simulink/SimMechanics based 
simulation that utilizes the Parameter Estimation Tool. 


Nomenclature 


C rotor = rotor damping constant 

CF = cost function for Newton’s least squares method 

/ m-rotor = moment of inertia for the MR (x-axis) rotor 

/ mt- rotor = moment of inertia for the MT (y-axis) rotor 

I rotor = generic rotor moment of inertia 

K rolgr = rotor spring constant 

M = generic data value (measured) 

S = generic data value (simulation) 

CO rotor = rotor natural frequency 

C rotor = rotor damping ratio 


I. Introduction 

Calculating an accurate nutation time constant (NTC), or nutation rate of growth, for a spinning upper stage is 
an important step for ensuring mission success. Spacecraft nutation, or wobble, is caused by energy dissipation 
anywhere in the system. Propellant slosh in the spacecraft fuel tanks is the primary source for this dissipation and, if 
it is in a state of resonance, the NTC can become short enough to violate mission constraints 1 . The Spinning Slosh 
Test Rig (SSTR), developed by NASA and Southwest Research Institute (SwRI), is a forced motion spin table 
where fluid dynamic effects in full-scale fuel tanks can be tested in order to obtain key parameters used to calculate 
the NTC. This is accomplished by independently varying nutation frequency vs. the spin rate and measuring force 
and torque responses on the tank. It has been used to predict parameters for the Genesis, CONTOUR, and STEREO 
missions whose tanks were mounted outboard from the spin axis. These parameters are incorporated into a 
mathematical model which uses mechanical analogs, such as pendulums and rotors, to simulate the force and torque 
resonances associated with fluid slosh. There is a high degree of uncertainty in predicting the effect of liquid 
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propellant motion in spinning spacecraft. Purely analytical methods of predicting the influence of onboard liquids 
have been generally unsatisfactory. The NTC values provided analytically are quite often significantly different than 
actual flight values. Hence, there is a need to identify conditions of resonance between the spacecraft nutation 
frequency and its liquid modes and to understand the general characteristics of the liquid motion. 

Extensive analysis has been done on the different tank shapes and locations, as well as the use of propellant 
management devices (PMD). A summary of these analyses, like that reported by Hubert 2 shows the vast differences 
in possible behaviors of different tank designs. For the off-spin-axis-mounted, cylindrical tanks with hemispherical 
end-caps that have been popular in a number of spacecraft programs, a number of relatively simple mechanical 
models have been developed. Hubert also notes that one of the most difficult aspects of employing such mechanical 
models is in the selection of appropriate parameters in the model. Use of mechanical analogs such as pendulums 
and rotors to simulate sloshing mass is a common alternative to fluid modeling. A homogeneous vortex model of 
liquid motion in spinning tanks and an equivalent mechanical rotor model was developed by Dodge 3 . An 
approximate theory of oscillations that predict the characteristics of the dominant inertial wave oscillation and the 
forces and moments on the tank are described. 

Most recently, the SSTR has been modified to simulate the on-axis spin motions of the centerline mounted tank 
used in the Pluto New Horizons (PNH) and Deep Impact (DI) spacecraft* as illustrated by the schematic diagram in 
Fig. 1. Diaphragm shapes, nutation frequencies, and fill levels were varied at a constant tank spin rate of =60.5 
RPM to examine how the various force and torque activity will influence the mechanical analog parameters. 



Figure 1. Schematic Diagram of On-Axis SSTR with DI Stand-off. 

Data is recorded while each fill level/diaphragm shape combination undergoes a nutation sweep where specific 
nutation frequencies range from 10-90% of the total spin rate. Two cameras were also added, along with particles in 
the liquid, to observe how the fluid moves within the tank for each nutation frequency. After testing, it was 
determined that the pendulum analog parameters could be set to zero because very little resonance activity was 
observed in the force response along the tank wall, leaving the rotor analog available to be the sole influencer of any 
torque resonance response. Resonance activity was observed in the torque response. This resonance activity, 
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caused by internal inertial wave oscillations, was verified by observing the particle motions in the tank with the on- 
board cameras. As a result of these observations, SwRI decided to model this resonance activity with two rotors. 
These ideal, or massless, rotors have inertia only about their primary axis of spin. Their associated parameters 
include inertia and spring/damping constants. The current method to identify these parameters involves laboriously 
hand-deriving equations of motion for both the SSTR and mechanical analog(s) and, by trial-and-error, comparing 
their results with experimental results. A strong desire exists to automate this method so that different analogs and 
parameter combinations can be quickly tested and results compared with measured data. The method presented 
accomplishes this task by using a modular SSTR MATLAB Simulink/SimMechanics-based simulation that 
identifies the parameters using the Parameter Estimation Tool 5 . 

A. Mechanical Analog Slosh Model Background 

Mechanical analogs, such as pendulums and rotors, are a common method of modeling fluid slosh. There are 
three main types of fluid slosh that can be observed when a PMD is not present (Fig. 2). These are surface wave, 
bulk fluid motion, and vortex mode slosh. Surface slosh and bulk fluid motion impose a force and torque response 
on the tank walls while the subsurface 
inertial waves only impose a torque response 
on the tank walls. In our test cases, a PMD 
(elastomeric diaphragm) was added and the 
tank rotation was restricted to being on-axis. 

The surface wave slosh became negligible 
and the bulk fluid motion was greatly 
reduced as indicated by a lack of resonance 
in the measured forces along the tank wall. 

This eliminated the need for a pendulum 
model to simulate force modes since none 
were present. However, resonance was 
taking place in the measured torque, 
accompanied by visible inertial waves. A 
rotor analog situated about the radial and 
tangential axis of the tank was developed to simulate these torque resonances. Figure 3 shows a diagram of one of 
these rotors. The rotor’s natural frequency and damping ratio can be calculated using the inertia and spring/damping 
constants as shown in Eqs. (1) and (2). 




(liquid reorients itself 
as tank rocks slowly) 



Vortex Mode 

(subsurface 
inertial waves ) 


Figure 2. Common Types of Liquid Slosh. 



Source: SwRI 

Figure 3. Schematic Diagram of Rotor Model. 
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An important factor that affects the fluid involves the 
shape of the diaphragm. Since the fill fraction of the tank 
varies from 30-80% during the testing, the diaphragm can 
take a wide variety of shapes. Bumps, ridges, troughs, 
and shifts in the diaphragm can have a significant effect 
on how the vortex mode interacts with the tank wall. 

During the tests, the diaphragm was manipulated into 
three standard shapes (mountain, crater, and yin-yang) as 
illustrated in Fig. 4. While there are many more Mountain Crater Yin -Yang 

diaphragm shapes possible, these three distinct shapes 

were felt to adequately encompass the overall shape Figure 4. Diaphragm Shapes Tested by SwRI. 
envelope. 

C. SwRI Data Analysis 

In the SSTR, the tank under test is suspended by three separate 3-axis force transducers. The data collected by 
these transducers ultimately gets resolved into two force and torque vs. time datasets about the X (radial) and Y 
(tangential) axis. The force about the X and Y axis is referred to as FX and FY while the torque is referred to as 
MR and MT respectively. The data supplied to the SSTR is the upper/lower motor RPM as well as the data 
collection start/stop time. Force, torque, upper/lower frame angular velocity, and upper/lower frame angle are 
recorded by the SSTR. All of the data is collected at 200Hz. 

Each test consists of a series of nutation sweeps where the nutation frequency (ratio of upper motor RPM to 
total RPM) varies from approximately 10-90% of the total spin rate. Each nutation sweep test usually consists of 9- 
14 individual nutation frequencies. Each test is designed to capture 30 nutation cycles with test times ranging from 
approximately 33 to 293 seconds. SwRI determined parameters for 16 nutation sweep tests as illustrated in Table 1. 
Once the tests are complete, a frequency domain analysis is performed and the signal amplitude and phase are 
plotted against a normalized nutation frequency. The rotor parameters are then systematically manipulated by hand 
until a best bit match with the measured data occurs. These parameters include the rotor inertia, spring stiffness 
coefficient, and damping coefficient. As Fig. 5 illustrates, there is no resonance activity in either the force or 
moment during a proof load test while there is resonance taking place in the radial torque during the tank test. The 
phase is defined as the difference between the force/torque output and the upper frame angular velocity. This 
relationship is used because the upper frame angular velocity is mechanically driven by the SSTR and does not 
change, regardless of the state of the liquid inside the tank. This provides the baseline for measuring the phase 
difference at any given nutation frequency. With this measurement methodology in place, the rotor parameters can 
now be adjusted to provide a best fit to the torque resonance as shown in Fig. 6. 



Table 1. PNH/DI On-Axis SSTR Test Matrix. 



Tank Volume Fill Fraction 

Mission 

Diaphragm Shape 

0.8 

0.706* 

0.624** 

0.5 

0.4 

PNH 

Crater 

X 

X 

X 

X 

X 

Mountain 

X 


X 

X 

X 

Yin- Yang 

X 


X 

X 

X 

Dl 

Crater 


X 




Mountain 


X 




Yin- Yang 


X 




* 

Deep Impact 

Flight F 

II Level 



Pluto New Horizons Fill Level 
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Figure 5. Processed Force, Torque, and Phase Data during Nutation Sweep. 
(62.4% Fill Level, Mountain Diaphragm Shape) 
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Figure 6. Liquid Analog Torque and Phase Responses during Nutation Sweep. 
(70.6% Fill Level, Crater Diaphragm Shape) 
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II. Modeling the SSTR Using MATLAB Simulink/SimMechanics 


In 2006, a lateral fuel slosh parameter estimation project to estimate pendulum parameters using the automated 
MATLAB Simulink/SimMechanics based parameter estimation method was completed 6 . Based on that project, it 
was decided to apply that modeling methodology to the SSTR and compare the results with the tried-and-tested 
results from SwRI. This method used the Simulink Parameter Estimation Tool. This tool provides a graphical 
interface where several optimization algorithms can be selected. The simulation incorporates the same slosh analogs 
used by SwRI, two torque rotors about the X (radial) and Y (tangential) axis respectively, so that direct comparisons 
between the hand-derived and automated methods can be made. 

A. Modeling the SSTR 

The Simulink/SimMechanics model represents the SSTR in such a way that the measured physical parameters 
can be directly translated to model parameters. These parameters include tank mass, moment of inertia(s), 
upper/lower frame spin RPM, physical geometry, CG locations, etc... The model was built using information 
supplied by SwRI and underwent an extensive calibration process while maintaining a high level of adjustability to 
accommodate future tests. Figure 7 gives an overall view of the model user interface. 



Figure 7. On-axis SSTR Simulation User Interface. 

The model can be divided up into three overall groups. These are, in order, the SSTR Drive Assembly, Spin 
Table and Tank Mount, and Fuel: Liquid Analog. Each of the first two groups defines the operating environment for 
the Fuel: Liquid Analog, and will remain the same regardless of the mechanical slosh analog used. This modular 
setup lets the user “drag and drop” many different fuel analogs into the model if desired. The SSTR Drive 
Assembly contains blocks that represent the upper and lower motors of the SSTR whose rotational speed can be 
varied independently from one another. Also, the initial starting position of the upper motor can be supplied, which 
is essential for matching to the physical setup. The tilt angle of 3.1 degrees is also defined here and provides the 
reference angle for the rest of the simulation. For example, if the SSTR were modified to have a different nutation 
angle, only one value would need to be changed in the simulation and the rest of the block orientations would reflect 
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the change automatically. One variable also controls whether or not the simulation is on-axis or off-axis. Currently 
it is set to zero for on-axis. The last variable that is included in the drive assembly is the distance to the spin center 
value. For PNH, it is set to 19.920 in. The drive assembly is setup to receive variable RPM vs. Time information if 
desired. 

The Spin Table and Tank Mount represent the second group and is the most complex of the three. The frame of 
the SSTR is accurately represented here and includes three sensors for the tank mount that are in the same location 
in space and orientation as their real-world counterparts. SimMechanics tends not to be stable when models have 
multiple supports since they are, technically, redundant. That is, if all three tank mounts were welded (zero degrees 
of freedom) to the upper portion of the frame, which is then welded to the lower portion of the frame, SimMechanics 
would only consider one weld sufficient and ignore the other two. To counteract this, special joints were added to 
support the virtual tank mount, causing it to float/rotate on/about each sensor’s local X-Y axis in order to record 
accurate data. The simulation is perfectly stable with this configuration. However, there is a very small transient 
region for about the first 0.08 seconds of the simulation. Since all of the estimation tests vary from roughly 33 to 
293 seconds, this transient can be ignored. 

The local coordinate systems for each of the three transducers are aligned so that their respective X-Y axis is 
tangent/perpendicular to an imaginary circle centered about them. This was done to prevent singularities from 
occurring in the revolute joints for each sensor. Since the SSTR sensors are all orientated the same way, two of the 
sensors’ FX and FY force datasets are converted in real-time using Simulink so that they align in common 
orientation identical to the SSTR. All of this happens automatically in the model and only requires interaction from 
the user if the sensor spacing is changed. For now, the sensor spacing is set at +/- 119.69 degrees*. Figure 8 
illustrates the overall layout of the model. 


Tank Mount Moment Sensor 



* Regardless of the sensor spacing, the forces will all add up to the same final value in the common orientation. 
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Unlike the force sensors, the moment sensor is located at the tank mount CG to provide the tangential and radial 
torque data, or MT and MR torque data respectively. Since SwRI has gone through great lengths to resolve the 
moments about this point, determining how to recombine the moment data from its three respective parts only to 
recombine it again for direct comparison to the measured data was deemed to be counter-productive. Handles to 
incorporate the force and moment DC offsets, or overall signal offsets, are also integrated into the spin table model 
if required. Currently, each of the three force sensors use one third of the value given by SwRI so that they all add 
up to a the final value while the moment sensor directly uses the given value. For the tests conducted, the DC bias 
was removed from the data in post processing so this feature was kept inactive. However, in keeping the model 
flexible for future use, the ability to incorporate DC bias was retained. 

The Fuel: Liquid Analog group is the simplest, yet most important part of the entire model. What makes this 
group unique is that the other two groups remain the same for each test (ignoring factors such as DC offsets and 
RPM values) while the liquid analog parameters will be actively changing during the parameter estimation. 
Modeling the rotor analogs used by SwRI in Simulink/SimMechanics is relatively straightforward as illustrated in 
Fig. 9. 



Figure 9. SwRI Rotor Analogs in Modeled in SimMechanics. 

Each massless § rotor is restricted to 1-DOF motion about the X (radial) and Y (tangential) axis and they only have 
inertia about their axis of spin as shown in Eq. 3. Attached to each rotor axle joint is a rotary spring/damper. If 


§ For model stability, each rotor was assigned an arbitrary mass of 0.01 lbm. 
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desired, each rotor can also start with user defined initial conditions. Each rotor joint is attached to the frozen 
liquid/tank frame. This contains the mass and moments of inertia of the fuel and the bare (empty) tank. 


Total Rotor Inertia 


I MR-Rotor 0 0 

0 I MT -Rotor 0 

0 0 0 


(3) 


The other component is the fixed, or non-sloshing, mass of the liquid along with the mass of the test tank and 
structure suspended within the force transducer array. Since the fuel mass was precisely measured by SwRI for each 
test, it was assumed to be known and was not included in the estimation. The calculated fixed mass moments of 
inertia for the proof-load tests were used as starting values and the products of inertia, being negligible, were set to 
zero. Table 2 shows the SSTR bare tank properties. For any given estimation, these are the minimum values for the 
fixed mass moments of inertia. 


Table 2. Bare SSTR Tank Mass Properties. 


Mass (Ibm) 

Ixx (lb*in 2 ) 

lyy (lb*in 2 ) 

Izz (lb*in 2 ) 

127.6 

15287 

12312 

22420 


In summary, six sets of measured data can be supplied to the simulation for parameter estimation. These 
include forces (FX/FY) along the tank wall and radial/tangential torques (MR/MT) resolved to the tank CG, as well 
as angular velocities of the tank CG (QX/QY) about each axis. Figure 10 compares the simulation results to the 
measured data when hard coding in the SSTR operating parameters provided by SwRI. The model’s FX, FY, MR, 
and MT outputs are a close match to the measured data and the angular velocities are almost identical to the 
measured data. This shows that the simulation output matches the SSTR measured output when provided the 
published operating parameters from SwRI for a given case. 
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Figure 10. Measured vs. Simulation using Values Provided by SwRI. 
(Mountain Diaphragm Shape at 62.4% Fill Level with a Nutation Frequency of 0.603 Hz) 
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B. Data Processing 

One of the fundamental differences between the SwRI and MATLAB methods is how the raw SSTR data is 
processed. Each sensor transmits raw data at 200Hz. SwRI immediately converts this data to the frequency domain 
where the amplitude, phase, and frequency are calculated to create the plots illustrated in Fig. 6. Since 
SimMechanics operates in the time domain, an alternative method had to be developed. First, the raw data is 
reduced to 50Hz in order to facilitate faster simulation times. Next, the data is run through a 5Hz low pass filter to 
remove noise and spikes. This is acceptable since the nutation frequency is less than 1Hz by default. As a final 
step, the DC bias is removed in order to reduce the number of parameters required to run the simulation. However, 
the handles to accommodate a DC bias were left in if needed. 

In order to conduct an accurate comparison between the two methods, Fourier tools were developed to convert 
the simulation data into the frequency domain. SwRI calculates this for 30 cycles at each nutation frequency, with 
time given at the beginning of each test to eliminate any transient region. A similar method is used for processing 
the simulation data. Observations show that the simulation transient region only lasts for the first few nutation 
cycles, so the method developed removes the first 10 nutation cycles to eliminate any transient region. This method 
is in agreement to the data processing done at SwRI. An added benefit of converting from the time domain is that it 
gives a better understanding of the signal components as illustrated in Fig. 11. With these tools in place, phase vs. 
nutation frequency/total spin rate plots can now be created. 


Frequency Spectrum Comparison Between Proof Load and Tank Data 



III. Tank Data Parameter Estimation with SwRI Rotors 


A. Simulink Parameter Estimation with Multiple Datasets 

The Parameter Estimation Tool provides a graphical interface that allows the user to choose from several 
optimization algorithms. In this case, the algorithm chosen is Newton’s method of non-linear least squares. This 
method modifies the model parameters and compares the resulting simulation data to the measured data in the time 
domain. The main objective pursued by the estimator is to lower the overall cost function, as defined in Eq. 4. A 
parameter, in this case, can be any variable in the MATLAB workspace. 

CF = Yj [abs(M ) - abs(S)] 2 (4) 

One of the more powerful features of the Parameter Estimation Tool is the ability to optimize parameters to 
multiple sets of measured data. Initially, the SSTR model was configured to optimize to measured data for one test. 
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Later, it was configured to accept data for tests covering the entire nutation sweep. Figure 12 illustrates the 
Parameter Estimation user interface showing an optimization termination notification for a run that covered seven 
sets of measured data (Test 003-Test 010). Tests 001 and 002 were not included due to the fact that the low nutation 
frequencies are outside the range of interest. Also, including these lengthy datasets would add a large amount of 
extraneous data that would needlessly increase the amount of time to perform each estimation. 



Figure 12. Parameter Estimation Tool Graphical User Interface. 

The estimator requires that the model must be prepared in a particular way to be able to run multiple dataset 
optimizations. First, the estimator requires that the desired output data be exported from the model using a Simulink 
“outport”. This data is then compared to the transient, or measured, data supplied to the estimator. Any 
combination of measured/transient data can be used as long as there is a corresponding outport block. The 
simulation of the SSTR currently contains six outport blocks as seen in Fig. 7. These are: Total FX, Total FY, Total 
MR, Total MT, Omega X, and Omega Y. These blocks can be connected into any part of the simulation to relay 
data back to the estimator. For example, if data from a single force transducer was the only set available, an outport 
block can be connected in the simulation and the estimator will use that data for the estimation. Several 
combinations ranging from using one dataset to using all six datasets were conducted to assess the best test 
combination for estimating the desired parameters. For the parameter estimation tests, Total FX, Total FY, Total 
MR, and Total MT datasets were used. The processed measured data is imported as a measured value vs. time 
dataset. Each dataset must start and end at the same time and have the same number of data points. 

Next, the parameters used by the SSTR Drive Assembly need to supplied to the estimator. This is accomplished 
by the use of “inports”. Just as with an outport, an inport can be located anywhere within a simulation. With the 
inports in place, the estimator can then supply the proper values required to produce the desired output data for each 
test. The data required to run the simulation are the lower frame RPM, upper frame RPM, and the upper frame 
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initial position. Each of these parameters is a scalar value supplied by SwRI for each test. If desired, the simulation 
can incorporate non-scalar input values for tests involving, for example, a non-constant upper/lower frame RPM. 

Once the measured data is entered, the parameters to be estimated can be defined. Any variable in the 
MATLAB workspace can be defined as a parameter. After an extensive investigation on the affects each parameter 
has in the resulting output data, a parameter estimation list was created as shown in Table 3. Large limit ranges 
were chosen for all but the rotor inertia parameters (in order to keep their values from becoming unreasonable). 
Only the rotor inertia values were restricted to a limiting case with the other parameters allowed an almost unlimited 
range of values to converge upon. It can be noted that common parameters can be left in the simulation as-is for 
each test because, as in the real SSTR, they remain static. These parameters include things like the tank structure 
mass, SSTR component mass properties, frame tilt angle, etc... 

Table 3. Parameters Used in Estimation Runs. 


Parameter Name 

Description 

Ixx 

Fixed Mass Ixx 

iyy 

Fixed Mass Iyy 

Izz 

Fixed Mass Izz 

rbxl 

MR Rotor Damping Constant 

rby2 

MT Rotor Damping Constant 

rkxl 

MR Rotor Spring Constant 

rky2 

MT Rotor Spring Constant 

ronelxx 

MR Rotor Inertia 

rtwolyy 

MT Rotor Inertia 


One of the benefits of the Parameter Estimation Tool is that the parameter estimation can be viewed in real- 
time. That is, the user can observe the parameters changing values as the optimization progresses and see how they 
affect the simulation response. This information is referred to by MATLAB as the measured vs. simulated response 
and the parameter trajectory. Figure 13 illustrates how the parameters change with each iteration while Fig. 14 
displays how the optimized solution compares with the measured data in a multiple dataset estimation. 



Figure 13. Example of a Parameter Trajectory for an Optimization. 


Once the parameter trajectories become constant or a convergence is observed visually, the simulation can be 
manually terminated by the user if desired and the parameter values from the last iteration will be retained. For 
more accurate estimations, the optimizer can iterate until a pre-defmed stopping criterion is met as indicated in the 
status window of the Parameter Estimator interface. If a more accurate solution is desired, the stopping criterion can 
be adjusted. In practice, however, it was found that increasing the stopping criterion settings greatly increased CPU 
time while only providing parameter values that lower the cost function by a small amount. 
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Simulation Data 
Measured Data 


Time (sec) 

Multiple Dataset Optimization: Before 

Measured vs Simulated Responses 



Multiple Dataset Optimization: After 


Figure 14. Before and After Comparison of Simulation and Measured Data. 

B. Test Summery and Results 

SwRI determined fixed mass inertia and rotor parameters for a total of 16 tank configurations. Each nutation 
sweep test consisted of 9-14 nutation frequency holds at frequencies ranging from roughly 10-90% of the total spin 
rate of -60.5 RPM. By estimating to the entire set of measured data, the rotor parameters can be optimized to 
reproduce the torque resonance response over the entire range of tested nutation frequencies. In order to verify the 
MATLAB approach, 3 cases were picked where noticeable resonance was observed as shown in Table 4. For each 
case, the estimator was given four sets of rotor limit ranges in order to determine if the estimator can converge closer 
to the measured data. The origins of this assumption are related to the particle observations in the fluid. During the 
spin tests, it was observed that all of the fluid in the tank was oscillating at certain nutation frequencies. It was 
decided to explore how this large vortex mode may influence the measured torques by modeling rotors with larger 
than average inertias. Each case first uses hard coded values for the rotor inertia range. The values used are those 
presented by SwRI. This range set is used as a benchmark to compare the two methods. In range sets two through 
four, the rotor inertia range is allowed to increase until it converges on a specific value. A fifth range set was added 
to the 70.6% Crater Diaphragm case that uses four measured data sets around a particular nutation frequency range 
using SwRI rotor inertias. This test was designed to observe how the rotor natural frequency can be manipulated by 
controlling the measured data processed and to observe any decrease in processing time needed to converge on a 
solution. The limits, and resulting rotor inertia values, are displayed in Table 5. 
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Table 4. Test Matrix for MATLAB Simulink/SimMechanics Method. 



Tank Volume Fill Fraction 

Mission 

Diaphragm Shape 

0.8 

0.706* 

0.624** 

0.5 

0.4 

PNH 

Crater 

X 

X 

X 

X 

X 

Mountain 

X 


X 

X 

X 

Yin-Yang 

X 


X 

X 

X 

Dl 

Crater 


X 




Mountain 


X 




Yin-Yang 


X 





Cases Chosen for MATLAB Method 


* 

Deep Impact Flight Fill Level 

** 

Pluto New Horizons Fill Level 


Table 5. Rotor Inertia Limit Summary. 



Limit Range 

Rotor Inertia (lb-in 2 ) 

Minimum 

Maximum 

Final Value 

Mountain 
Diaphragm at 
62.4% Fill Level 

SwRI 

MR 

- 

- 

700 

MT 

- 

- 

700 

Limit #1 

MR 

100 

1000 

1000 

MT 

100 

1000 

1000 

Limit #2 

MR 

100 

1200 

1200 

MT 

100 

1200 

1200 

Converged 

MR 

100 

5000 

1272.5 

MT 

100 

5000 

1895.7 

Crater Diaphragm 
at 70.6% Fill Level 

SwRI 

MR 

~ 

- 

600 

MT 

- 

~ 

500 

Limit #1 

MR 

100 

1000 

1000 

MT 

100 

1000 

1000 

Limit #2 

MR 

100 

1500 

1500 

MT 

100 

1500 

1500 

Converged 

MR 

100 

15000 

7007.4 

MT 

100 

15000 

8024.3 

SwRI (T6-T9) 

MR 

- 

- 

600 

MT 

- 

~ 

500 

Mountain 
Diaphragm at 
80.0% Fill Level 

SwRI 

MR 


- 

1500 

MT 

- 

- 

1400 

Limit #1 

MR 

500 

2000 

2000 

MT 

500 

2000 

2000 

Limit #2 

MR 

500 

2500 

2500 

MT 

500 

2500 

2500 

Converged 

MR 

500 

20000 

6969.4 

MT 

500 

20000 

7361 .4 


As the table illustrates, the rotor inertia consistently hits the upper limits during all but the last optimization. 
Correspondingly, the fixed mass inertia decreases as the rotor inertia increases to maintain conservation of 
momentum as shown in Table 6. Table 7 shows how the cost function decreases as the rotor inertia constraints are 
increased. 
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Table 6. Summary of Liquid Fixed Mass Inertias. 



Limit Range 

Ixx* (lb-in 2 ) 

lyy* db-in 2 ) 

Izz* (lb-in 2 ) 

Mountain Diaphragm 
at 62.4% Fill Level 

SwRI 

3825.9 

4313.1 

4665.8 

Limit #1 

3615.6 

4102.2 

4767.4 

Limit #2 

3464.3 

3950.2 

4823.6 

Converged 

3329.3 

3083.7 

4777.4 

Crater Diaphragm at 
70.6% Fill Level 

SwRI 

6647.1 

7226.7 

6740.1 

Limit #1 

5665.5 

6153.0 

6453.9 

Limit #2 

4742.2 

5259.5 

6254.0 

Converged 

572.8 

0.0 

7588.4 

SwRI (T6-T9) 

6827.6 

7478.8 

7153.3 

Mountain Diaphragm 
at 80.0% Fill Level 

SwRI 

4116.9 

4775.9 

4013.4 

Limit #1 

3360.0 

3895.2 

3967.7 

Limit #2 

2635.1 

3177.3 

3923.5 

Converged 

0.0 

125.6 

5873.0 

* In addition to the bare SSTR tank inertia 



Table 7. Cost Function Summary. 



Limit Range 

Cost Function 

Mountain Diaphragm 
at 62.4% Fill Level 

SwRI 

388370 

Limit #1 

386510 

Limit #2 

385730 

Converged 

380230 

Crater Diaphragm at 
70.6% Fill Level 

SwRI 

670630 

Limit #1 

604510 

Limit #2 

555650 

Converged 

424960 

SwRI (T6-T9) 

268450 

Mountain Diaphragm 
at 80.0% Fill Level 

SwRI 

881770 

Limit #1 

818100 

Limit #2 

773240 

Converged 

643410 


The rotor spring and damping value limits are set to be large enough to accommodate any rotor inertia value. 
This was done to allow the rotor natural frequency to maintain a specific natural frequency, regardless of the rotor 
inertia. Table 8 contains the results for the rotor spring and damping constants for each test. As the table shows, the 
rotor natural frequencies are very close to those calculated by SwRI. The two exceptions are the 70.6% case (which 
uses all of the measured data) and the cases where the rotors were allowed to converge on a maximum value. The 
cause for the rotor natural frequencies starting out low in the 70.6% case lies in the nature of how the optimizer 
converges on a set of measured data. Since each test is 30 nutation cycles, the length of the low nutation frequency 
tests is much longer. This means that there are more data points in the lower nutation frequency region, potentially 
causing the optimizer to bias results to the left. Since the phase of the low nutation frequency measured data is at 
least half that of the peak resonance frequency, the lowest cost function happens to occur when the phase resonance 
of the simulation is at a lower frequency. As the rotor inertia increases, the simulation can better match all of the 
datasets, moving the rotor peak offset back to the right. These shifts in rotor natural frequency can be mitigated by 
choosing data around the nutation frequency region where it is expected that resonance takes place as the results for 
the 70.6% T6-T9 case show. An added benefit for ignoring data outside of the nutation frequency region of interest 
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is that it speeds up the estimation since there is less data to process. For the cases where the rotor inertia upper limit 
cap was lifted, the natural frequency tended to increase to a higher frequency. For the Rotor Inertia Limit #1 and #2 
tests, the rotor natural frequency does increase by a small amount, but not nearly as much as for the converged rotor 
cases. For all of the tests, the spring constant increases with rotor inertia to maintain the natural frequency as 
governed by Eq. 1. While the present difference between the damping constants is large for some of the tests, it must 
be noted that the values for damping are very small. 

Table 8. Summary of Rotor Spring and Damping Constants. 



Spring Constant 
(ft.lb/deg) 

Damping Constant 
(ft.lb/deg/sec) 

Natural F 
_lh 

requency 

z) 

Mountain Diaphragm at 62.4% Fill Level 

MR 

MT 

MR 

MT 

MR 

MT 

SwRI Results 

0.0519 

0.0519 

0.0026 

0.0026 

0.7066 

0.7066 

MATLAB Results: SwRI Rotor Inertias 

0.0550 

0.0543 

0.0024 

0.0024 

0.7266 

0.7224 

Percent Difference from SwRI Results 

5.84% 

4.61% 

6.00% 

7.37% 

2.84% 

2.24% 

MATLAB Results: Rotor Inertia Limit #1 

0.0823 

0.0807 

0.0049 

0.0048 

0.7440 

0.7366 

MATLAB Results: Rotor Inertia Limit #2 

0.1021 

0.0994 

0.0069 

0.0068 

0.7562 

0.7463 

MATLAB Results: Converged Rotor Inertias 

0.1095 

0.1658 

0.0094 

0.0196 

0.7604 

0.7668 


Crater Diaphragm at 70.6% Fill Level 

MR 

MT 

MR 

MT 

MR 

MT 

SwRI Results 

0.0328 

0.0274 

0.0019 

0.0016 

0.6066 

0.6066 

MATLAB Results: SwRI Rotor Inertias 

0.0186 

0.0151 

0.0024 

0.0015 

0.4569 

0.4504 

Percent Difference from SwRI Results 

43.22% 

44.83% 

24.65% 

5.47% 

24.69% 

25.75% 

MATLAB Results: SwRI Rotor Inertias (T6-T9) 

0.0328 

0.0274 

0.0013 

0.0009 

0.6066 

0.6070 

Percent Difference from SwRI Results (T6-T9) 

0.11% 

0.19% 

29.78% 

41.97% 

0.00% 

0.06% 

MATLAB Results: Rotor Inertia Limit #1 

0.0325 

0.0314 

0.0053 

0.0051 

0.4677 

0.4597 

MATLAB Results: Rotor Inertia Limit #2 

0.0501 

0.0482 

0.0104 

0.0101 

0.4740 

0.4649 

MATLAB Results: Converged Rotor Inertias 

0.4467 

0.5224 

0.1337 

0.1710 

0.6546 

0.6616 


Mountain Diaphragm at 80.0% Fill Level 

MR 

MT 

MR 

MT 

MR 

MT 

SwRI Results 

0.0766 

0.0765 

0.0046 

0.0044 

0.5858 

0.6061 

MATLAB Results: SwRI Rotor Inertias 

0.0784 

0.0708 

0.0053 

0.0044 

0.5926 

0.5831 

Percent Difference from SwRI Results 

2.35% 

7.39% 

15.85% 

0.26% 

1.15% 

3.79% 

MATLAB Results: Rotor Inertia Limit #1 

0.1046 

0.1034 

0.0086 

0.0084 

0.5929 

0.5895 

MATLAB Results: Rotor Inertia Limit #2 

0.1314 

0.1274 

0.0128 

0.0127 

0.5943 

0.5854 

MATLAB Results: Converged Rotor Inertias 

0.4573 

0.4863 

0.0760 

0.0881 

0.6642 

0.6664 


Once the estimations were complete, plots similar to Figs. 5 and 6 were created to see how the model reacts to 
the different parameters. Figures 15 and 16 show the force amplitude over a nutation frequency sweep for the 
62.4% Mountain Diaphragm case. Since the torque rotors have no effect on the forces, the amplitude remains the 
same regardless of the rotor size or state. The force amplitude plots are similar for the other two cases. Figures 1 7 
through 22 show the moment amplitude over a nutation frequency sweep for each case. While the overall 
magnitude of the torques during each rotor limit test are the same for the 62.4% case, the torque magnitude using 
converged rotor values grows to be well outside of family in the 70.6% and 80.0% cases. 

Figures 23 through 28 show the phase offset plots for each case. As the plots show, the model data comes 
closer to matching the measured data as the rotor inertias are allowed to increase. The only dataset that does not 
match well is the phase plot for the SwRI rotor inertias in Figures 25 and 26 for the 70.6% case due to the lower 
rotor natural frequencies. However, the phase offset can be moved to the correct location as demonstrated using the 
T6-T9 case. In addition to limiting to datasets around the desired nutation frequency, the rotor natural frequency can 
be manually restricted by the user. However, for the purposes of these tests, it was left entirely to the estimator to 
choose the rotor parameters that would yield the lowest cost function. 
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Figure 15. FX Amplitude During Nutation Frequency Sweep (62.4% Mountain). 
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Figure 16. FY Amplitude During Nutation Frequency Sweep (62.4% Mountain). 
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Figure 17. MR Amplitude During Nutation Frequency Sweep (62.4% Mountain). 
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Figure 18. MT Amplitude During Nutation Frequency Sweep (62.4% Mountain). 
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Figure 19. MR Amplitude During Nutation Frequency Sweep (70.6% Crater). 
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Figure 20. MT Amplitude During Nutation Frequency Sweep (70.6% Crater). 
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Figure 21. MR Amplitude During Nutation Frequency Sweep (80.0% Mountain). 
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Figure 22. MT Amplitude During Nutation Frequency Sweep (80.0% Mountain). 
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Figure 23. MR to ftX Phase During Nutation Frequency Sweep (62.4% Mountain). 



- *♦ - Measured 

SwRI Rotor Inertias 

Rotor Inertia Limit #1 

Rotor Inertia Limit #2 

Converged Rotor Inertias 



Figure 24. MT to 1XY During Nutation Frequency Sweep (62.4% Mountain). 
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Figure 25. MR to ftX During Nutation Frequency Sweep (70.6% Crater). 
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Figure 26. MT to ftY During Nutation Frequency Sweep (70.6% Crater). 
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Figure 27. MR to 12X During Nutation Frequency Sweep (80.0% Mountain). 
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Figure 28. MT to DY During Nutation Frequency Sweep (80.0% Mountain). 
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IV. Discussion of Results and Current Work 


After estimating parameters for these three cases, it can be seen that the MATLAB Simulink/SimMechanics 
model does an adequate job at matching to the measured data. This is significant, considering the entirely different 
approach from the SwRI method and the parameter selection is 100% automated. With the rotor inertias hard coded 
to match the SwRI values, the spring and damping constants converged to values close to those supplied by SwRI in 
the 62.4% and 80.0% Mountain Diaphragm cases. For the 70.6% Crater Diaphragm case, the parameters adequately 
matched after optimizing to data only around the nutation frequencies of interest. An added benefit of using data 
only near the nutation frequencies of interest is that the Parameter Estimation Tool can run faster since it no longer 
has to process extraneous data. As with the SwRI method, the estimator tended to converge on similar rotor spring 
and damping constants in order to recreate the proper phase offsets throughout each nutation frequency sweep. 
Early on in the testing it was found that the estimator tended to increase the rotor inertia and take away from the 
fixed mass inertia, resulting in a lower cost function. As the phase plots show, the simulation data does match the 
phase offset of the measured data better as the rotor inertias increase. This suggests that most of the fluid is rotating, 
creating a large subsurface inertial wave. Physically, this is what appears to be happening in the suspended particle 
observation videos from the SSTR on-board cameras. It is logical to assume that larger rotors may be able to more 
accurately simulate the measured torque response up to a certain point. The amplitude of the torque data must be 
taken into account when increasing the rotor inertias because it may deviate as it did for the 70.6% and 80.0% cases 
using the converged rotor parameters. This large torque amplitude deviation suggests that there is an upper limit to 
the rotor inertias, even though the cost function is lower and the phase offset is a much better match than the other 
cases. For the Limit #1 and #2 cases, when the rotor inertias were allowed to increase by a small amount, the phase 
offsets better matched the data with only a small increase in the torque amplitude. The trade-offs between 
optimizing to torque amplitude vs. optimizing to phase offset deserve further study. 

The MATLAB method of estimating slosh analog parameters has proved to be a good tool for augmenting 
SwRI’s method. There are several benefits to using a Simulink/SimMechanics based SSTR model. Most notable, 
liquid slosh analogs can be quickly constructed in the model to investigate how the forces and torques react as 
different mechanical analogs are tried out. This eliminates the need to tediously hand-derive the equations of 
motion for each analog during the initial development phase. Once a final slosh analog is chosen, the model can 
then be used to validate the hand-derived equations. The test rig portion of the model also behaves very much like 
the real SSTR so possible modifications and/or upgrades can be tested before any actual hardware is replaced. 
Finally, the simulation inherently uses the entire signal frequency spectra to identify paramaters. In the future, this 
capability could be used to identify slosh parameters for several modes at once. 

The slosh analog can also be easily moved to Simulink-based spacecraft models to directly calculate the NTC. 
Efforts are currently underway to create a PNH NTC model to incorporate the SSTR Fuel: Liquid Analog block as 
illustrated in Figure 29. An NTC estimate will be generated by this model using the various rotor parameters. 
These results will then be compared to flight data. This technique was used for the NTC models of the Genesis, 
CONTOUR and STEREO spacecraft respectively. The match of predicted NTC to flight NTC got better with each 
subsequent use. 



Figure 29. Simulink/SimMechanics-Based PNH NTC Model. 
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Figure 30. Liquid Analog Torque and Phase Responses during Nutation Sweep. 
(62.4% Fill Level, Mountain Diaphragm Shape) 
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Figure 31. Liquid Analog Torque and Phase Responses during Nutation Sweep. 
(80.0% Fill Level, Mountain Diaphragm Shape) 
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